clear all
set more off

global resultspath /INSERTPATH

local stackedgraphs
local incomes high middle low 
local temperatures 10 20 30

foreach income of local incomes {
foreach temperature of local temperatures {

if "`income'" == "high" {
local income_level 10.819
}

if "`income'" == "middle" {
local income_level 9.21
}

if "`income'" == "low" {
local income_level 6.907
}

local countrysetname 	all	
local lhsvar		ln_revenue_worker 	/*ln_revenue ln_revenue_worker
					numberofemployees ln_numberofemployees*/
local sector 		manu		  /*all manu serv*/
local transform 	poly4	   	/*bins poly4 polyabove1 polyabove2 poly2*/
local int			tripleinteracted	  /*interacted uninteracted tripleinteracted*/
local ln_gdppc		`income_level'
local avg_tmax		`temperature'
local CI		yes		/*yes no*/
local starttemp		-10
local endtemp		40
local time_fe 		country_year_fe /*state_year_fe*/
local tosave		yes		/*yes no*/
local capital		no
local weather_controls			/*_precip_only 	_precip_plus_lags   _temp_lags*/

clear 
local omit 15
set obs 61
gen temperature = 0.5
keep temperature
replace temperature = -15+_n-1
if "`regressiontype'" == "bins" {
replace temperature = temperature+0.5
}

estimate use "$resultspath/ster/`countrysetname'_`sector'_`transform'_`int'_`lhsvar'_`weights'_`time_fe'_`weather_controls'.ster"

gen ln_revenue_worker = .
gen numberofemployees = .
gen ln_numberofemployees = .

if "`regressiontype'" == "poly2" & "`int'" == "interacted" {
predictnl yhat = temperature*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax') ///
+temperature^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax') ///
-(`omit'*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax') ///
+`omit'^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax')) ///
, se(se) ci(lowerci upperci)
}
if "`regressiontype'" == "poly3" & "`int'" == "interacted" {
predictnl yhat = temperature*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax') ///
+temperature^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax') ///
+temperature^3*(_b[poly_3]+_b[c.poly_3#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_3#c.avg_tmax_lr]*`avg_tmax') ///
-(`omit'*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax') ///
+`omit'^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax') ///
+`omit'^3*(_b[poly_3]+_b[c.poly_3#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_3#c.avg_tmax_lr]*`avg_tmax')) ///
, se(se) ci(lowerci upperci)
}
if "`regressiontype'" == "poly4" & "`int'" == "interacted" {
predictnl yhat = temperature*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax') ///
+temperature^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax') ///
+temperature^3*(_b[poly_3]+_b[c.poly_3#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_3#c.avg_tmax_lr]*`avg_tmax') ///
+temperature^4*(_b[poly_4]+_b[c.poly_4#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_4#c.avg_tmax_lr]*`avg_tmax') ///
-(`omit'*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax') ///
+`omit'^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax') ///
+`omit'^3*(_b[poly_3]+_b[c.poly_3#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_3#c.avg_tmax_lr]*`avg_tmax') ///
+`omit'^4*(_b[poly_4]+_b[c.poly_4#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_4#c.avg_tmax_lr]*`avg_tmax')) ///
, se(se) ci(lowerci upperci)
}

if "`regressiontype'" == "poly4" & "`int'" == "tripleinteracted" {
predictnl yhat = temperature*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_1#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax') ///
+temperature^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_2#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax') ///
+temperature^3*(_b[poly_3]+_b[c.poly_3#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_3#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_3#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax') ///
+temperature^4*(_b[poly_4]+_b[c.poly_4#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_4#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_4#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax') ///
-(`omit'*(_b[poly_1]+_b[c.poly_1#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_1#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_1#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax') ///
+`omit'^2*(_b[poly_2]+_b[c.poly_2#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_2#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_2#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax') ///
+`omit'^3*(_b[poly_3]+_b[c.poly_3#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_3#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_3#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax') ///
+`omit'^4*(_b[poly_4]+_b[c.poly_4#c.ln_gdp_pc]*`ln_gdppc'+_b[c.poly_4#c.avg_tmax_lr]*`avg_tmax' ///
+_b[c.poly_4#c.ln_gdp_pc#c.avg_tmax_lr]*`ln_gdppc'*`avg_tmax')) ///
, se(se) ci(lowerci upperci)
}

if "`regressiontype'" == "poly4" & "`int'" == "uninteracted" {
predictnl yhat = temperature*_b[poly_1]+temperature^2*_b[poly_2]+temperature^3*_b[poly_3]+temperature^4*_b[poly_4] ///
-(`omit'*_b[poly_1]+`omit'^2*_b[poly_2]+`omit'^3*_b[poly_3]+`omit'^4*_b[poly_4]), se(se) ci(lowerci upperci)
}

if "`regressiontype'" == "bins" & "`int'" == "interacted" {
predictnl yhat = (_b[bin_5C_35_Inf]+_b[c.bin_5C_35_Inf#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_35_Inf#c.avg_tmax_lr]*`avg_tmax')*(temperature==37.5) ///
+(_b[bin_5C_30_35]+_b[c.bin_5C_30_35#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_30_35#c.avg_tmax_lr]*`avg_tmax')*(temperature==32.5) ///
+(_b[bin_5C_25_30]+_b[c.bin_5C_25_30#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_25_30#c.avg_tmax_lr]*`avg_tmax')*(temperature==27.5) ///
+(_b[bin_5C_20_25]+_b[c.bin_5C_20_25#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_20_25#c.avg_tmax_lr]*`avg_tmax')*(temperature==22.5) ///
+0*(temperature==17.5) ///
+(_b[bin_5C_10_15]+_b[c.bin_5C_10_15#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_10_15#c.avg_tmax_lr]*`avg_tmax')*(temperature==12.5) ///
+(_b[bin_5C_5_10]+_b[c.bin_5C_5_10#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_5_10#c.avg_tmax_lr]*`avg_tmax')*(temperature==7.5) ///
+(_b[bin_5C_0_5]+_b[c.bin_5C_0_5#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_0_5#c.avg_tmax_lr]*`avg_tmax')*(temperature==2.5) ///
+(_b[bin_5C_n5_0]+_b[c.bin_5C_n5_0#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_n5_0#c.avg_tmax_lr]*`avg_tmax')*(temperature==-2.5) ///
+(_b[bin_5C_n10_n5]+_b[c.bin_5C_n10_n5#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_n10_n5#c.avg_tmax_lr]*`avg_tmax')*(temperature==-7.5) ///
+(_b[bin_5C_nInf_n10]+_b[c.bin_5C_nInf_n10#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.bin_5C_nInf_n10#c.avg_tmax_lr]*`avg_tmax')*(temperature==-12.5) ///
, se(se) ci(lowerci upperci)
drop if yhat == 0 & temperature != 17.5
}

if "`regressiontype'" == "bins" & "`int'" == "uninteracted" {
predictnl yhat = _b[bin_5C_35_Inf]*(temperature==37.5) ///
+_b[bin_5C_30_35]*(temperature==32.5) ///
+_b[bin_5C_25_30]*(temperature==27.5) ///
+_b[bin_5C_20_25]*(temperature==22.5) ///
+0*(temperature==17.5) ///
+_b[bin_5C_10_15]*(temperature==12.5) ///
+_b[bin_5C_5_10]*(temperature==7.5) ///
+_b[bin_5C_0_5]*(temperature==2.5) ///
+_b[bin_5C_n5_0]*(temperature==-2.5) /// 
+_b[bin_5C_n10_n5]*(temperature==-7.5) ///
+_b[bin_5C_nInf_n10]*(temperature==-12.5) ///
, se(se) ci(lowerci upperci)
drop if yhat == 0 & temperature != 17.5
}

if "`regressiontype'" == "polyabove1" & "`int'" == "interacted" {
predictnl yhat = ((temperature-30)*(_b[polyAbove30_1]+_b[c.polyAbove30_1#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.polyAbove30_1#c.avg_tmax_lr]*`avg_tmax'))*(temperature>30) ///
+((5-temperature)*(_b[c.polyBelow5_1]+_b[c.polyBelow5_1#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.polyBelow5_1#c.avg_tmax_lr]*`avg_tmax'))*(temperature<5) ///
, se(se) ci(lowerci upperci)
}

if "`regressiontype'" == "polyabove1" & "`int'" == "tripleinteracted" {
predictnl yhat = ((temperature-30)*(_b[polyAbove30_1]+_b[c.polyAbove30_1#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.polyAbove30_1#c.avg_tmax_lr]*`avg_tmax' + _b[c.polyAbove30_1#c.avg_tmax_lr#c.ln_gdp_pc]*`avg_tmax'*`ln_gdppc'))*(temperature>30) ///
+((5-temperature)*(_b[c.polyBelow5_1]+_b[c.polyBelow5_1#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.polyBelow5_1#c.avg_tmax_lr]*`avg_tmax' + _b[c.polyBelow5_1#c.avg_tmax_lr#c.ln_gdp_pc]*`avg_tmax'*`ln_gdppc'))*(temperature<5) ///
, se(se) ci(lowerci upperci)
}

if "`regressiontype'" == "polyabove1" & "`int'" == "uninteracted" {
predictnl yhat = ((temperature-30)*_b[polyAbove30_1])*(temperature>30) ///
+((5-temperature)*_b[polyBelow5_1])*(temperature<5), se(se) ci(lowerci upperci)
}

if "`regressiontype'" == "polyabove2" & "`int'" == "interacted" {
predictnl yhat = ((temperature-30)^2*(_b[polyAbove30_2]+_b[c.polyAbove30_2#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.polyAbove30_2#c.avg_tmax_lr]*`avg_tmax'))*(temperature>30) ///
+((5-temperature)^2*(_b[c.polyBelow5_2]+_b[c.polyBelow5_2#c.ln_gdp_pc]*`ln_gdppc' ///
+_b[c.polyBelow5_2#c.avg_tmax_lr]*`avg_tmax'))*(temperature<5) ///
, se(se) ci(lowerci upperci)
}

if "`regressiontype'" == "polyabove2" & "`int'" == "uninteracted" {
predictnl yhat = ((temperature-30)^2*_b[polyAbove30_2])*(temperature>30) ///
+((5-temperature)^2*_b[polyBelow5_2])*(temperature<5), se(se) ci(lowerci upperci)
}

if "`lhsvar'" == "ln_revenue" {
local varname revenue
}
if "`lhsvar'" == "ln_revenue_worker" {
local varname rvw
}
if "`lhsvar'" == "numberofemployees" {
local varname numemp
}
if "`lhsvar'" == "ln_numberofemployees" {
local varname ln_emp
}

if "`lhsvar'" == "ln_revenue_worker" {
graph tw rarea upperci lowerci temperature if temperature > `starttemp', color(ltblue) || ///
line yhat temperature if temperature > `starttemp', lcolor(black) ///
yscale(range(-0.01 0.005)) ylabel(-0.01(0.005).005, nogrid) xscale(range(`starttemp' `endtemp')) ///
xlabel(`starttemp'(10)`endtemp') graphregion(fcolor(white)) xtitle("Daily Maximum Temperature") ytitle("% Revenue Per Worker") ///
name("lnGDP_`income'_temp_`temperature'_`varname'") legend(off)
}

if "`lhsvar'" == "ln_revenue" {
graph tw rarea upperci lowerci temperature if temperature > `starttemp', color(ltblue) || ///
line yhat temperature if temperature > `starttemp', lcolor(black) ///
yscale(range(-0.015 0.005)) ylabel(-0.015(0.01)0.005, nogrid) xscale(range(`starttemp' `endtemp')) ///
xlabel(`starttemp'(10)`endtemp') graphregion(fcolor(white)) xtitle("Daily Maximum Temperature") ytitle("% Revenue") ///
name("lnGDP_`income'_temp_`temperature'_`varname'") legend(off)
}

if "`lhsvar'" == "ln_numberofemployees" {
graph tw rarea upperci lowerci temperature if temperature > `starttemp', color(ltblue) || ///
line yhat temperature if temperature > `starttemp', lcolor(black) ///
yscale(range(-0.005 0.005)) ylabel(-0.005(0.005)0.005, nogrid) xscale(range(`starttemp' `endtemp')) ///
xlabel(`starttemp'(10)`endtemp') graphregion(fcolor(white)) xtitle("Daily Maximum Temperature") ytitle("% {&Delta} Employment") ///
name("lnGDP_`income'_temp_`temperature'_`varname'") legend(off)
}

if "`lhsvar'" == "numberofemployees" {
graph tw rarea upperci lowerci temperature if temperature > `starttemp', color(ltblue) || ///
line yhat temperature if temperature > `starttemp', lcolor(black) ///
yscale(range(-6 2)) ylabel(-6(2)2, nogrid) xscale(range(`starttemp' `endtemp')) ///
xlabel(`starttemp'(10)`endtemp') graphregion(fcolor(white)) xtitle("Daily Maximum Temperature") ytitle("{&Delta} Employment") ///
name("lnGDP_`income'_temp_`temperature'_`varname'") legend(off)
}

local stackedgraphs `stackedgraphs' lnGDP_`income'_temp_`temperature'_`varname'
}
}

graph combine `stackedgraphs' 

if "`tosave'" == "yes" {
graph export "$resultspath/Graphs/`countrysetname'_`sector'_`transform'_`int'_`lhsvar'_`weights'_`time_fe'_`weather_controls'.png", replace width(16000)
}



